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I. INTRODUCTION. 



Black hole simulations are known to provide a severe test of Numerical Relativity, even 
in the one-dimensional (ID) case [Q]. These strong field scenarios imply a wide dynamical 
range involving very different time and length scales. The coordinate degrees of freedom 
must be used with special care in order to prevent the numerical code to crash when space- 
time approaches a singularity. This is the idea underlying the different singularity avoiding 
conditions which have been proposed and tested in Numerical Relativity [0-^. This ap- 
proach, however, is not free from problems. Spacetime dynamics gets locally frozen near the 
singularity, leading to an abrupt transition zone around the horizon. This can cause either 
steep space gradients |^ or spikes in the radial metric function p| which make difficult to 
maintain the accuracy or even lead to code crashing during numerical evolution. 

A major step towards a singularity-proof scheme in Numerical Relativity is the use of 
a horizon boundary condition. The idea |^|^ is to reduce drastically the dynamical range 
by evolving just the observable region and imposing a suitable boundary condition on or 
slightly inside the horizon, which is a one-way membrane. This idea has been actually 
implemented in Ref. ^ by the combined use of a "horizon locking" coordinate system and a 
"causal" finite difference scheme which is very similar to the "causal reconnection" scheme 
introduced in Ref. Numerical evolution was there shown to proceed without significant 
errors beyond the limit oft = 100m, where the previous codes that used fixed boundaries |l[ 
became unstable or extremely inaccurate. The results presented in Ref. [|I| correspond to the 
maximal slicing condition 0. This does not mean that that code could not use alternative 
slicings, like the harmonic one. It is simply due to the fact that maximal slicing happens 
to be very robust, leading to longer numerical black hole evolution than other singularity 
avoiding slicings 

Singularity avoidance, however, is not our only guideline in constructing a Numerical 
Relativity code. We want to use a system of evolution equations which actually translates 
the causal structure of the spacetime: it should be a hyperbolic system of partial differential 
equations with the local speed of light as characteristic speed. It has been shown recently 
|]T^,|TT| that the use of a harmonic time coordinate (harmonic slicing [Q) leads to such 
hyperbolic evolution systems, which can also be expressed in flux-conservative form. In 
that way, we can get the well known structure of the hydrodynamic equations and the 
huge arsenal of the Computational Fluid Dynamics (CFD) methods is at our disposal. The 
straightforward use of such methods has yet resulted into a 3D numerical code which is 
able by now to evolve vacuum spacetimes admitting periodic boundary conditions. 

In the present work, we will use CFD methods to improve the quality of harmonic slicing 
black hole codes. The flrst part (sections 2-4) of the paper deals with the standard approach, 
where we use a flnite difference discretization in a flxed grid. The horizon boundary condition 
is implemented in the second part (sections 5-8) by using a moving numerical grid. The use 



of a moving grid is not new in Numerical Relativity. Wilson [T^ used it to deal with 
the Relativistic Hydrodynamics equations. The same "adaptive mesh" technique was then 
extended to the fleld equations and successfully applied to the study of axisymmetric stellar 
collapse 0,|l5|. Here we provide a new application of this technique to the study of black 
hole evolution, where we use the grid speed to keep the numerical mesh outside the horizon. 
This approach is compared with the "horizon locking" condition introduced in Ref. M. 
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II. THE ID BLACK HOLE. 



Let us write the spherically symmetric (ID) vacuum line element in isotropic coordinates: 

which is locally isometric to the Schwarzschild metric. The lapse function 

a=P^^ (2) 

vanishes at p = m/2, which corresponds to the black hole horizon (r = 2m in standard 
Schwarzschild coordinates). 

The form of the line element is convenient for numerical applications: (i) it is easy 
to express ([^) in Cartesian coordinates; this allows one to consider the same problem as a 
test for generic (3D) numerical codes, and (ii) the metric (|lD is invariant under the inversion 
symmetry 

P . m/2 



m/2 p 

mapping the black hole outer region (p > m/2) into the inner part and vice versa; this pro- 
vides a consistent inner boundary condition at the inversion point p = m/2 (the "throat"), 
allowing one to evolve just the exterior part of the black hole spacetime ("wormhole" evo- 
lution: see Ref. |jl|). 

We shall use the space part of the line element (|1|) to provide the initial data for our 
numerical models. The singular lapse function (j^) will be however replaced by a non-singular 
symmetric function with initial value 

a{0,p) = constant (4) 

in order to get a different spacetime slicing so that the metric components are no longer static 
but evolving in time. The line element is expressed in the generic spherically symmetric form 

rfs2 _ _^2^^^ ^^2 ^ ^2^^^ ^^2 ^ Y\t, p) dVt. (5) 

The initial condition @) preserves the inversion symmetry (^ of the spacetime. This 
means that the throat connecting the two isometric sheets will remain fixed at p = m/2 and 
inner boundary conditions can be imposed there consistently at every step of the evolution 
process. Note however that the black hole horizon |TB[ moves with the local light speed 



c = a/X, (6) 

which will be now different from zero everywhere. This means that the horizon will start 
moving outwards, away from the fixed throat position. The coincidence between throat and 
horizon in the static form (|l|) was because the singular lapse choice (§) caused the vanishing 
of the local speed of light at that point. 
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In order to locate the black hole horizon [0, we will make use of the apparent horizon 
local definition 



Yja + Y'jX = 0, (7) 

where dots and primes stand for time and space derivatives, respectively. The mass function 
M{t,p), defined by 

2M/Y = 1 + [{Y/aY - iY'/X)% (8) 

gives at every instant the total amount of mass enclosed by a sphere of coordinate radius p. 
Of course, it does coincide with the constant parameter m in our case. Condition (|^ then 
implies that Y = 2m at the black hole horizon. 



III. THE NUMERICAL CODE. 

Our code is a finite difference version of Einstein field equations in first order form. This 
means that the set of basic quantities includes not only the lapse function a and the spatial 
metric components (the contravariant ones g'^^ in our case), but also their first order space 
and time derivatives. The time derivative of the lapse function is given by imposing the time 
coordinate to be harmonic (harmonic slicing [^]). For the sake of simplicity, the remaining 
set of independent quantities will be described in what follows as the components of a single 
vector valued function u. These quantities are taken to be independent because we consider 
the constraint equations as first integrals of our evolution system, which are imposed on the 
initial data only (free evolution approach). The conservation of the constraints can then be 
used as an accuracy test. 

Under these conditions, it has been shown [|10| that in the generic 3D case the evolution 
system can be written as an hyperbolic system of balance laws 

dtU + dkF\u) = S{u), (9) 

where the fiuxes F'^ and sources S are vector valued functions of u. Moreover, it has been 
shown [0] that there is an infinite family of hyperbolic evolution systems which share only 
the physical solutions (the ones actually satisfying the constraint equations). We will use 
here a spherically symmetric (ID) version of one of these systems, which is explicitly given 
in the Appendix A. 

The source terms in (|^) take into account the nonlinear part Einstein's equations. We 
have used an operator splitting approach by considering separately the source driven evolu- 
tion 



and the transport process 



dtu = S{u) (10) 



dtu + dpF{u) = 0. (11) 



These different processes are then combined (Strang splitting |T^) to obtain a second order 
accurate discretization of the full equation (Mj. 
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In order to analyze the importance of the source terms in the overall evolution, we 
have implemented both a standard second order Runge-Kutta along with a sophisticated 
high order Bulirsch-Stoer method |T^. We have found that the accurate modeling of the 



transport step was far more significant than the source treatment in the overall evolution. 

We have used a second order upwind method (see Appendix B) to deal with the transport 
step. Let us note that the hyperbolicity of the system ([TT|) allows one to express it as two 
uncoupled subsystems with the structure of (the first order form of) the wave equation |]TU 



This is useful when implementing the upwind method because one can easily construct the 
linear combinations of the original variables which propagate along light rays going in the 
forward or backward directions (characteristic variables, see the Appendix). 

The boundary conditions at the inner boundary are then to be imposed on the forward 
propagating combinations only. This is done by allowing for the inversion symmetry (|^) at 
the throat (we will give more details in Sec. 5). Conversely, the outer boundary conditions 
will affect only to the backward propagating combinations. 



IV. FIXED GRID RESULTS. 

We have performed our computations with an evenly spaced grid of 200 points ranging 
from 1 to 40 Schwarzschild radii, thus obtaining a resolution that can be already tested in 
3-D codes thus allowing future comparison of results. The accuracy was monitored by 



computing the mass function (H), which will keep its constant value M{t, p) = m only if the 
constraint equations are preserved. We have found that this provides an extremely sensitive 
error test. 

The results presented in Fig. |l| actually show the mass function, as computed from the 
numerically evolved quantities (m = 2 in our case). The time values correspond to the 
proper time of the outermost evolving point. Errors are big around the horizon position, 
which is computed from (^ and marked with a cross, in spite of the huge dynamical range 
one gets at the black hole throat as one approaches the singularity, as it is clearly shown in 
Fig. ^ for the radial metric component. 

This apparent paradox is explained by the collapse of the speed of light in the inner 
zone which locally freezes the evolution there, as it is clearly shown in Fig. ^ This behavior 
amounts to the well known collapse of the lapse which is generic to singularity avoiding 
coordinate systems, like harmonic or maximal slicing. It allows us to deal with the huge 
dynamical range in the metric components near the singularity (see Fig. but it has a 
perverse side effect: the appearance of steep gradients near the horizon which are actually 
the main source of numerical errors. 

One can of course try to reduce these errors by increasing either the grid resolution or 
the accuracy of the numerical method. We have rather preferred, in keeping with Ref. to 
avoid the gradients by evolving just the exterior part of the black hole spacetime (a moving 
boundary problem) with the moving grid techniques sketched in the following section. 
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V. A MOVING GRID APPROACH. 



Let us introduce a new radial function r, which is related to the previous one p in a time 
dependent way: 

P = f{t,r), (12) 

so that a grid of fixed points with respect to the new variable r will be moving with respect 
to the original one, attached to p. The first order derivatives of (12) 

V = dj, A = drf (13) 

can be interpreted as the relative speed and dilation factor, respectively, as computed by an 
observer attached to the moving grid. 

The structure of the set of conservation laws (p!TD does not change when written in terms 
of the new radial function 



dt (A u) + dr [F{u) -Vu] = 0, (14) 

where we have transformed the quantities it in a scalar way. The new characteristic matrix 
keeps the same set of eigenvectors as the original one, so that hyperbolicity is preserved. 
The new characteristic speeds are obtained from the old ones simply by subtracting the grid 
speed V and then dividing by the dilation factor A. Notice that the change of variables (|T2p 
introduces just one new degree of freedom, the dilation factor being obviously related to the 
grid speed: 

dtA - drV = 0. (15) 

Even the static choice V = gives us the possibility of choosing a nontrivial dilation 
factor to deal with numerical grids which are not evenly spaced with respect to the original 
variable p. We actually used it to obtain the fixed grid results of the previous section. The 
reason has to do with the treatment of the inner boundary: the image under (|^) of an evenly 
spaced set of points is no longer evenly spaced. Equations (p!4D with V = provided an 
elegant way to preserve the overall accuracy at the inner boundary [T^ without running into 
stability problems. 

The next simplest case is the linear one, where the dilation factor depends on the time 
coordinate only. This leads to a linear speed profile which can be determined by the boundary 
values. One can apply this to the black hole problem by demanding the grid speed to coincide 
with the local speed of light (characteristic speed) at the horizon. In that way, one is placing 
the inner boundary at the horizon (instead of at the throat), which will remain attached to 
a specific node of the moving grid. The profile can be fully determined by demanding the 
grid speed to vanish at (or near to) the external boundary. We have found this characteristic 
boundary approach very convenient to avoid the steep gradient zone behind the horizon wile 
keeping a high degree of accuracy in the black hole exterior, as we show in the following 
section. 
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VI. MOVING GRID RESULTS. 



We have redone with the moving grid the same computations presented in Sec. 4, under 
the same initial conditions. The hnear relationship between the fixed and moving grid 
coordinates allows us to plot our moving grid results in terms of the original variable p, 
allowing a direct comparison with Sec. 4. In that sense, Fig. | is to be compared with the 
previous Fig. |l|. The successive plots start now at the horizon position because the inner 
part is no longer computed. Note that the error in the mass has decreased drastically: this 
confirms that the larger errors in Fig. |I] were due to the steep gradient zone behind the 
horizon. 

Note that we are now placing 200 evenly spaced points in the region between the horizon 
and the outer boundary. The evolution can then be pursued until this region gets very 
small, leading to an extremely low dilation factor (extremely high characteristic speed) 
which freezes the evolution. In this moving grid approach, the lifetime of a numerical black 
hole is just the time it takes the horizon to arrive to the grid outer boundary ||2^ . 

We have also tested other initial values for the lapse function, different from (^). We 
found that the fixed grid code crashed in some cases, due to large errors, while the moving 
grid one kept its low error profile in every case. We conclude that the moving grid approach 
is both a more robust and more accurate way to deal with the ID black hole problem, 
significantly improving the accuracy of the results in every case. 



VII. RELATED APPROACHES. 

A similar approach, which would be more appealing to relativists, is to consider the 
transformation (0) as a change of spacetime coordinates, not just affecting the grid motion, 
but also the metric quantities. The line element (^) would then transform into: 

ds'^ = -a^{t, r) dt^ + X2(t, r) {^dr + V dtf + Y'^{t, r) dU^ (16) 

so that a radial shift vector appears: 

13' = V/A (17) 

and the radial metric component gets multiplied by the dilation factor. 

The new degree of freedom can now be used in a different way, by imposing a fixed form 
for the radial metric function and adding the shift (3 to the list of the variables The 
constant value of the radial metric function ensures that the use of a non-uniform dilation 
factor does not lead to a premature freezing of the evolution. But the main price to pay is 
that the structure of the original system (|Tl|) is modified at the risk of losing hyperbolicity 
and, with that, the possibility of applying standard numerical algorithms. We are currently 
studying the structure of the general 3+1 (or ADM) evolution system in order to see whether 
or not hyperbolicity can actually be preserved when introducing shift vectors and/or when 
using other algebraic gauge conditions. 

There are of course other adaptive grid approaches, based on mesh refinement algorithms 



21|, which have been successfully applied to the ID scalar field case |22]. Such powerful 
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methods include sophisticated bookkeeping routines to manage the grid, putting more reso- 
lution just where it is needed, without modifying the equations in any way. The application 
of these methods to the generic 3D case is a big challenge to the present day computing and 
programming capabilities. 

VIII. DISCUSSION AND OUTLOOK. 

We are aware that many of the difficulties one would encounter in evolving really dynamic 
spacetimes are avoided in the Schwarzschild test. As a first example, even in the ID case, 
let us remember that a black hole can be formed in the collapse of a supermassive star. The 
horizon forms at a given instant in the evolution and new horizons can appear later. The 
evolution of the inner boundary will be then piecewise continuous, jumping every time that 
a new outermost horizon appears. 

As stated in section 2, our code contains an horizon finding routine, based on the apparent 
horizon local definition (|^), which can detect when and where a new horizon does appear. 
This routine can also detect whether the horizon is actually ahead of the inner boundary due 
to cumulative numerical errors in computing the light speed there. We have also implemented 
another routine which, when the difference exceeds two grid zones, makes the inner boundary 
to jump from its previous location (the first grid point) to the actual horizon position, 
discarding the grid points between the old and new positions and recomputing the speed 
profile accordingly. We have found that a few of these jumps do not compromise the stability 
of the code and therefore we believe that the discontinuous horizon motion in ID dynamical 
spacetimes can be dealt with by using the same technique. 

One would encounter other difficulties when trying to extend this formalism to the mul- 
tidimensional black hole case. For instance, the horizon can expand in an anisotropic way 
so that the resulting speed profile can lead to a highly distorted numerical mesh. To solve 
this problem, one can make use of the relationship (0) between the grid speed and the shift 
vector. One could demand that the velocity field satisfies either the "minimal strain" or 
"minimal distortion" or a similar condition. The corresponding set of elliptic equations 
should be solved by using the light speed components as inner boundary values. 

But possibly the main difficulty in going to the 3D case is that the horizon finding routine 
can not be a straightforward generalization of the ID one. The multidimensional analogous 
of contains explicitly the unit normal to the horizon surface (something one does not 
know a priori). This is the heel of Achilles of 3D adaptive grid black hole codes. We recently 
heard about promising progress on that subject [^,|^, and it seems that it will help in 
solving such difficulty in a near future. 
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APPENDIX A: EVOLUTION EQUATIONS. 



The set of metric variables to evolve is 

u={c,rr,g'',g'',q\,q%,D\,D%) 

where C = ot^f^ is the local speed of light and we have introduced the shortcuts D'^ ^ = 
drlnig"^"^), D% = drln{g^^), Vr = D% - \D\ - U, and U = drln{a). 

With this notation, the vacuum Einstein Evolution Equations in spherical symmetry can 
be written in first order form as follows: 



ar. = c 



dtC = -c\% , 



2q\Lr + {q\ - q\)D\ 



dtg'' = Cg'W , 



dt {q\ + q)-dr [C{D\ + 2r, - 2L,)] = C {q\f - {D\f 



(Ala) 
(Alb) 
(Ale) 
(Aid) 
(Ale) 



dt {q\ + q)-dr [C{D\ - 2L,) 



C 



q\{2q\ - q\) + 2^ - {D'e f + 2L,D% 



, (Alf) 



dt{D\-2Lr)^dr[C{q\ + q)] , 



dt{D\-2Lr)^dr[C{q% + q) 



(Alg) 



(Alh) 



where we have also noted q = q^i = q^j. + 2q^g. 
The nontrivial characteristic quantities are 

w"^ = C(D\ + 2r, - 2Lr) ± {q\ + q) , < = C{D% - 2L,) ± {q% + q) , (A2) 

with characteristic speeds ±C. 



APPENDIX B: DISCRETIZATION OF THE TRANSPORT PROCESS. 

The finite difference version of the transport equation 

dtu + dpF{u) = 

is constructed in two steps. 



(Bl) 
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In the first step, one computes all the variables at the cell interfaces (points i ± |) at an 
intermediate time level + |. The resulting predictions 

are then used to compute the Fluxes 

Ftlli = F{<:ilh (B3) 

so that, in the second step, the next time level is computed in an explicit flux-conservative 
way: 

= < - i - Ftiji) ■ m 

Note that the stability of this explicit scheme implies that the time interval of every step 
is limited by the causality condition ensuring that the characteristic speeds are lower than 
one grid zone per timestep. 

Allowing for the fact that this second step is centered both in space and in time, first 
order accuracy in the interface values leads to second order accuracy in the final result, 
because the leading error terms cancel when subtracting the Fluxes between the right and 
left interfaces. 

In order to obtain the interface values, we first compute the forward and backward first 
order predictions for every quantity, 

„,FW 3..n l^.n I dt ( zpn rpn \ ('Dr^\ 

%-l/2 - 2 - 2^i-2 - 2d^ [■t i_i - -t i_2) , (Bo) 



Ui-i/2 - - 2%+i -2d^[Fi^^-Fi). (B6j 

Note that, when dealing with shocks or large gradients, these one-sided predictions are 
known to produce spurious oscillations: the predicted values at the i — \ interface may get 
out of the interval defined by the values at the grid points i and i + 1. This is anomalous 
in the sense that we are modeling a transport process and the causality condition ensures 
that nothing coming from outside of this interval can reach the interface at the intermediate 
time level. We detect this spurious behaviour when one of the anomalous sequences {uf, 
"^r+i' i'^f+i/2^ '^r^ "^r+i) monotomc. In these cases, we take iifj!f/2 — '^i+i 

'^f+1/2 ~ ''^^1 respectively, to make sure that our predictions do not get out of range. 

We perform then a local transformation at every interface from our original quantities 
u to the set characteristic variables w which are given in the preceding Appendix. The 
hyperbolicity of our transport equations ensures that this is a one to one invertible transfor- 
mation. The interface values for each characteristic component w are taken to be either the 
forward or backward prediction, depending on the sign of the corresponding characteristic 
speed (positive or negative, respectively). The interface values of the original quantities u 
are finally recovered by inverting the transformation. 
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FIGURES 



FIG. 1. Plots of the mass function (fixed grid case), which should be equal to 2 everywhere. 

The maximum error is around 15%. The successive horizon positions are marked with a cross. 
Note that the different plots coincide in the inner part, where evolution is freezing due to the light 
speed collapse (see Fig. 3). 

FIG. 2. Sucessive plots of the time evolution of the logarithm of the metric component (jrr- 
The crosses show the successive positions of the black hole horizon. Steep gradients appear just 
behind the horizon. The resulting dynamical range rises up to eighty orders of magnitude. 

FIG. 3. Plots of the coordinate light speed, showing a sudden collapse of the inner part (the 
black hole interior), due to the singularity avoidance of the harmonic slicing. 

FIG. 4. Same as Fig. 1 for the moving grid case. The maximum error is now less than 0.5%. 
Note that the scale has been enlarged by one order of magnitude. 
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